# Convenient jupyter setup
%load_ext autoreload
%autoreload 2
#%load_ext rpy2.ipython
import re
import pathlib
import geopandas as gpd
import pandas as pd
import rioxarray as rxr
from rasterio.enums import Resampling
# Visualisation imports
from holoviews import opts
import holoviews as hv
from holoviews.operation.datashader import rasterize
from holoviews.element import tiles
hv.extension("bokeh")
import numpy as np
from src.constants import (PAISAGENSLIDAR_PATH, EBALIDAR_PATH, DATA_PATH,
WGS84, WEBMERCATOR, SIRGAS_BRAZIL)
from src.utils.os import list_content, latest_matching_file
from src.visualisation.image_loader import fetch_metrics, visualise, load_image, load_raster, load_bounds, VISUALISATION_CACHE, stylise, create_raster_image
from src.visualisation.visualisation_params import *
# Global parameter whether to save
save = False
esri = tiles.ESRI().redim(x="Easting", y="Northing")
from src.data.jrc_loading import get_jrc_data
from src.processing.canopy_height_models import longo2016_acd, asner2014_acd, asner2014_colombia_acd
paisagens_master = gpd.read_feather(PAISAGENSLIDAR_PATH / "paisagens_master_table.feather")
eba_master = gpd.read_feather(EBALIDAR_PATH / "eba_master_table.feather")
#paisagens_master.iloc[0]
paisagens_survey_name = "SAN_A02_2014"
eba_survey_name = "NP_T-0859"
eba_survey_name2 = "NP_T-0408"
p_sample = paisagens_master.loc[paisagens_survey_name]
e_sample = eba_master.loc[eba_survey_name]
e_sample2 = eba_master.loc[eba_survey_name2]
samples = [p_sample, e_sample, e_sample2]
lon_min = min([sample.lon_min for sample in samples])
lat_min = min([sample.lat_min for sample in samples])
lon_max = max([sample.lon_max for sample in samples])
lat_max = min([sample.lat_max for sample in samples])
# Get bounds
l,b,r,t = gpd.GeoSeries(p_sample.geometry_bbox, crs=SIRGAS_BRAZIL).to_crs(WEBMERCATOR).bounds.values[0]
#xs, ys = np.array(list(gpd.GeoSeries(sample.geometry_metrics, crs=SIRGAS_BRAZIL).to_crs(WEBMERCATOR)[0].exterior.coords)).T
#poly = hv.Polygons({"x": xs, "y": ys}).redim(x="Easting", y="Northing").opts(alpha=0.2)
resolution = 50
chms = [load_raster(sample.name, "canopy_height_model", 1, from_crs=sample.crs, to_crs=None) for sample in samples]
chm_means = [chm.rio.reproject(chm.rio.crs, resolution=resolution, resampling=Resampling.average) for chm in chms]
# Filter high outliers
chm_means = [chm.where(chm < 100) for chm in chm_means]
longo2016_acds = [longo2016_acd(chm_mean.where(chm_mean < 100)) for chm_mean in chm_means]
asner2014_acds = [asner2014_acd(chm_mean.where(chm_mean < 100)) for chm_mean in chm_means]
asner2014_colombia_acds = [asner2014_colombia_acd(chm_mean.where(chm_mean < 100)) for chm_mean in chm_means]
transition_map = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="TransitionMap_MainClasses").squeeze().rio.reproject(WEBMERCATOR)
first_deforestation = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="DeforestationYear").squeeze().rio.reproject(WEBMERCATOR)
first_degradation = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="DegradationYear").squeeze().rio.reproject(WEBMERCATOR)
last_deforestation = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="LastDeforestationYear").squeeze().rio.reproject(WEBMERCATOR)
last_degradation = get_jrc_data(lon_min, lat_min, lon_max, lat_max, dataset="LastDegradationYear").squeeze().rio.reproject(WEBMERCATOR)
(visualise("SAN_A02_2014", "TransitionMap", 30, transition_map)
+ visualise(p_sample.name, "canopy_height_model", 1, datashader_aggregator="mean", from_crs=p_sample.crs)
+ visualise("SAN_A02_2014", "LastDeforestation", 30, last_deforestation.where(last_deforestation > 1) )
+ visualise("SAN_A02_2014", "FirstDeforestation", 30, first_deforestation.where(first_deforestation > 1))
+ visualise("SAN_A02_2014", "LastDegradation", 30, last_degradation.where(last_degradation > 1) )
+ visualise("SAN_A02_2014", "FirstDegradation", 30, first_degradation.where(first_degradation > 1))
).cols(1)
#get_jrc_data(*sample[["lon_min", "lat_min", "lon_max", "lat_max"]]).rio.reproject(WEBMERCATOR)
#deforestations = [deforestation.rio.reproject_match(chm_mean) for chm_mean in chm_means]
#degradations = [degradation.rio.reproject_match(chm_mean) for chm_mean in chm_means]
(do I need to account for disturbance as well? Check JRC user guide)
#years_since_degradation = [(sample.year - degradation.where(degradation > 0)) for sample, deforestation in zip(samples, deforestations)]
#years_since_deforestation = [(sample.year - deforestation.where(deforestation > 0)) for sample, deforestation in zip(samples, deforestations)]
#years_since_deforestation = [year.where(year > 0) for year in years_since_deforestation]
chm_means[0].rio.rep
CRS.from_epsg(32721)
import xarray as xr
from tqdm.autonotebook import tqdm
<ipython-input-46-1d2b4d910b30>:2: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console) from tqdm.autonotebook import tqdm
survey_name = "SAN_A02_2014"
survey = paisagens_master.loc[survey_name]
desired_resolution = 50
relevant_datasets = ("TransitionMap_MainClasses", "DeforestationYear", "DegradationYear",
"LastDeforestationYear", "LastDegradationYear")
# Load canopy height model (CHM)
chm = load_raster(survey.name, "canopy_height_model", 1, from_crs=survey.crs, to_crs=None)
# Aggregate CHM to desired resolution
chm_mean = chm.rio.reproject(chm.rio.crs, resolution=desired_resolution, resampling=Resampling.average)
# Extract relevant JRC datasets
jrc_data = {dataset: get_jrc_data(survey.lon_min,
survey.lat_min,
survey.lon_max,
survey.lat_max,
dataset=dataset
).squeeze()#.rio.reproject(survey.crs)
for dataset in relevant_datasets}
annual_change = xr.concat([get_jrc_data(survey.lon_min,
survey.lat_min,
survey.lon_max,
survey.lat_max,
dataset="AnnualChange",
year = year
).squeeze().assign_coords({"year": year})
for year in tqdm(range(1992, survey.survey_year + 1))],
dim="year")
# Get bounds
l,b,r,t = gpd.GeoSeries(survey.geometry_bbox, crs=SIRGAS_BRAZIL).to_crs(WEBMERCATOR).bounds.values[0]
xs, ys = np.array(list(gpd.GeoSeries(survey.geometry_metrics, crs=SIRGAS_BRAZIL).to_crs(WEBMERCATOR)[0].exterior.coords)).T
poly = hv.Polygons({"x": xs, "y": ys}).redim(x="Easting", y="Northing").opts(alpha=0.5)
from src.visualisation.visualisation_params import *
from bokeh.models import FixedTicker
from src.processing.jrc_processing import last_of_value
def last_observation(annual_change, jrc_class_value, first_observation=None):
if first_observation is not None:
last_observation = xr.zeros_like(first_observation)
else:
last_observation = xr.zeros_like(annual_change[0])
# Compute last year of observation (Note: Due to the size of the datasets this takes
# ~10 min, even with numba speedup)
last_observation.data = last_of_value(
annual_change.data, value=jrc_class_value
)
# Fix all negative/fill values and all values prior to 1990 with the first year of
# deforestation.
if first_observation is not None:
last_observation = xr.concat(
[first_observation, last_observation + 1990], dim="observations"
).max(dim="observations")
else:
last_observation = xr.concat(
[xr.zeros_like(last_observation), last_observation + 1990], dim="observations"
).max(dim="observations")
return last_observation.astype(float)
last_deforested = last_observation(annual_change, 3).rio.reproject(WEBMERCATOR)
last_deforested
<xarray.DataArray (y: 122, x: 232)>
array([[ 0., 0., 0., ..., 2012., 2012., 2000.],
[ 0., 0., 0., ..., 2012., 2012., 2012.],
[ 0., 0., 0., ..., 2012., 2012., 2012.],
...,
[2012., 2012., 2012., ..., 0., 2012., 2012.],
[2012., 2012., 2012., ..., 2012., 2012., 2012.],
[2012., 2012., 2012., ..., 2012., 2012., 2012.]])
Coordinates:
* x (x) float64 -6.113e+06 -6.113e+06 ... -6.106e+06 -6.106e+06
* y (y) float64 -3.652e+05 -3.652e+05 ... -3.688e+05 -3.688e+05
year int64 1992
band int64 1
spatial_ref int64 0
Attributes:
_FillValue: -9999.0
grid_mapping: spatial_refarray([[ 0., 0., 0., ..., 2012., 2012., 2000.],
[ 0., 0., 0., ..., 2012., 2012., 2012.],
[ 0., 0., 0., ..., 2012., 2012., 2012.],
...,
[2012., 2012., 2012., ..., 0., 2012., 2012.],
[2012., 2012., 2012., ..., 2012., 2012., 2012.],
[2012., 2012., 2012., ..., 2012., 2012., 2012.]])array([-6113084.994615, -6113054.983846, -6113024.973077, ..., -6106212.528475,
-6106182.517706, -6106152.506937])array([-365194.35894 , -365224.369709, -365254.380478, -365284.391247,
-365314.402016, -365344.412785, -365374.423555, -365404.434324,
-365434.445093, -365464.455862, -365494.466631, -365524.477401,
-365554.48817 , -365584.498939, -365614.509708, -365644.520477,
-365674.531246, -365704.542016, -365734.552785, -365764.563554,
-365794.574323, -365824.585092, -365854.595861, -365884.606631,
-365914.6174 , -365944.628169, -365974.638938, -366004.649707,
-366034.660476, -366064.671246, -366094.682015, -366124.692784,
-366154.703553, -366184.714322, -366214.725091, -366244.735861,
-366274.74663 , -366304.757399, -366334.768168, -366364.778937,
-366394.789706, -366424.800476, -366454.811245, -366484.822014,
-366514.832783, -366544.843552, -366574.854322, -366604.865091,
-366634.87586 , -366664.886629, -366694.897398, -366724.908167,
-366754.918937, -366784.929706, -366814.940475, -366844.951244,
-366874.962013, -366904.972782, -366934.983552, -366964.994321,
-366995.00509 , -367025.015859, -367055.026628, -367085.037397,
-367115.048167, -367145.058936, -367175.069705, -367205.080474,
-367235.091243, -367265.102012, -367295.112782, -367325.123551,
-367355.13432 , -367385.145089, -367415.155858, -367445.166627,
-367475.177397, -367505.188166, -367535.198935, -367565.209704,
-367595.220473, -367625.231243, -367655.242012, -367685.252781,
-367715.26355 , -367745.274319, -367775.285088, -367805.295858,
-367835.306627, -367865.317396, -367895.328165, -367925.338934,
-367955.349703, -367985.360473, -368015.371242, -368045.382011,
-368075.39278 , -368105.403549, -368135.414318, -368165.425088,
-368195.435857, -368225.446626, -368255.457395, -368285.468164,
-368315.478933, -368345.489703, -368375.500472, -368405.511241,
-368435.52201 , -368465.532779, -368495.543548, -368525.554318,
-368555.565087, -368585.575856, -368615.586625, -368645.597394,
-368675.608164, -368705.618933, -368735.629702, -368765.640471,
-368795.65124 , -368825.662009])array(1992)
array(1)
array(0)
mpl.cm.get_cmap('viridis', 30)
<matplotlib.colors.ListedColormap at 0x7ff00938e340>
viridisBig = mpl.cm.get_cmap('viridis', 512)
newcmp = ListedColormap(viridisBig(np.linspace(0.25, 0.75, 256)))
plot_examples([viridis, newcmp])
<matplotlib.colors.ListedColormap at 0x7ff01019c070>
(
hv.Image(annual_change.loc[survey.survey_year].rio.reproject(WEBMERCATOR)).opts(tools=["hover"],
**opt_annual_change,
**opt_size) * poly
+
hv.Image(annual_change.where(annual_change==4).loc[survey.survey_year].rio.reproject(WEBMERCATOR)).opts(tools=["hover"],
**opt_annual_change,
**opt_size) * poly
+
hv.Image(last_deforested.where(last_deforested > 0), vdims=["deforestation_year"]).opts(**opt_size,
colorbar=True,
cmap= mpl.cm.get_cmap('viridis', 30),
tools=["hover"]) * poly
).cols(1)
FixedTicker()
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-179-82ebe2fe7f42> in <module> ----> 1 FixedTicker({key: val["name"] for key, val in annual_change_legend.items()}) TypeError: __init__() takes 1 positional argument but 2 were given
def pixel_based_recovery(survey_name)
last_deforestation.rio.reproject_match()
(30.005401194024046, -30.005401194024046)
last_deforestation.rio.clip_box()
was_deforested = last_deforestation > 0
is_recovering
%%time
years_since_deforestation = []
for observation, chm_mean in zip(samples, chm_means):
print(observation.survey_year)
# Selectors
deforested = deforestation > 0
never_degraded = degradation == 0
not_distrubed_since_deforested = deforestation > degradation
deforested_before_observation = deforestation < observation.survey_year
# Convert deforestation year into recovery period
recovery_period = (observation.survey_year
- deforestation.where(deforested
& deforested_before_observation
& (never_degraded | not_distrubed_since_deforested)
)
)
# Reproject and match recovery period area to area with lidar observation
recovery_period = recovery_period.rio.reproject_match(chm_mean)
years_since_deforestation.append(recovery_period)
2014 2017 2016 CPU times: user 421 ms, sys: 20 ms, total: 441 ms Wall time: 447 ms
accumulated_heights = [chm_mean.where(recovery_period > 0) for chm_mean, recovery_period in zip(chm_means, years_since_deforestation)]
accumulated_biomasses_longo2016 = [biomass.where(recovery_period > 0) for biomass, recovery_period in zip(longo2016_acds, years_since_deforestation)]
accumulated_biomasses_asner2014 = [biomass.where(recovery_period > 0) for biomass, recovery_period in zip(asner2014_acds, years_since_deforestation)]
accumulated_biomasses_asner2014c = [biomass.where(recovery_period > 0) for biomass, recovery_period in zip(asner2014_colombia_acds, years_since_deforestation)]
xs = np.concatenate([recovery_period.squeeze().data.flatten() for recovery_period in years_since_deforestation])
ys_height = np.concatenate([height.squeeze().data.flatten() for height in accumulated_heights])
ys_longo2016 = np.concatenate([biomass.squeeze().data.flatten() for biomass in accumulated_biomasses_longo2016])
ys_asner2014 = np.concatenate([biomass.squeeze().data.flatten() for biomass in accumulated_biomasses_asner2014])
ys_asner2014c = np.concatenate([biomass.squeeze().data.flatten() for biomass in accumulated_biomasses_asner2014c])
# Visualisation imports
from holoviews import opts
import holoviews as hv
from holoviews.operation.datashader import rasterize
from holoviews.element import tiles
hv.extension("bokeh")
df = pd.DataFrame(np.vstack((xs,
ys_height,
ys_longo2016,
ys_asner2014,
ys_asner2014c)
).T,
columns=["recovery_period [years]",
"mean_height [m]",
"longo2016_biomass [kgC/m2]",
"asner2014_biomass [kgC/m2]",
"asner2014c_biomass [kgC/m2]"]
).dropna().reset_index(drop=True)
table = hv.Table(df)
points = hv.Points(df)
df.head()
| recovery_period [years] | mean_height [m] | longo2016_biomass [kgC/m2] | asner2014_biomass [kgC/m2] | asner2014c_biomass [kgC/m2] | |
|---|---|---|---|---|---|
| 0 | 12.0 | 15.705148 | 6.877455 | 9.425895 | 6.906647 |
| 1 | 12.0 | 14.114010 | 5.698736 | 8.514473 | 6.031612 |
| 2 | 17.0 | 17.577612 | 8.385388 | 10.492824 | 7.967184 |
| 3 | 17.0 | 14.914895 | 6.280082 | 8.973813 | 6.468921 |
| 4 | 12.0 | 11.979083 | 4.269915 | 7.283663 | 4.898956 |
hv.BoxWhisker(table, kdims=["recovery_period [years]"], vdims="mean_height [m]").sort().opts(**opt_size)
BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('index', 51), ('mean_height [m]', 51), ('mean_height_left_square_bracket_m_right_square_bracket', 0)
max_observation = int(points["recovery_period [years]"].max())
hist = points.hist(dimension=["recovery_period [years]"],
num_bins=max_observation,
bin_range=(0.5, max_observation + 0.5),
adjoin=False).opts(**opt_size, tools=["hover"])
hist
hv.BoxWhisker(table, kdims=["recovery_period [years]"], vdims="mean_height [m]").sort().opts(**opt_size)
BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('index', 20), ('mean_height [m]', 20), ('mean_height_left_square_bracket_m_right_square_bracket', 0)
points = hv.Points(df)
max_observation = int(points["recovery_period [years]"].max())
hist = points.hist(dimension=["recovery_period [years]"],
num_bins=max_observation,
bin_range=(0.5, max_observation + 0.5),
adjoin=False).opts(**opt_size, tools=["hover"])
hist
hv.BoxWhisker(table, kdims=["recovery_period [years]"], vdims="longo2016_biomass [kgC/m2]").sort().opts(**opt_size)
BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('index', 58), ('longo2016_biomass [kgC/m2]', 58), ('longo2016_biomass_left_square_bracket_kgC_over_m2_right_square_bracket', 0)
rasterize(points).opts(cmap="kbc_r", cnorm="eq_hist").opts(**opt_size, tools=["hover"])
In the JRC TMF data specification it says that the DeforestationYear marks the first year of deforestation. What we
need for our regrowth analysis, however, is the last year of deforestation. This will tell us from when on the forest was undistrubed and could resume its growth.
JRC procedure: The JRC-TMF classification process starts out by mapping disturbances in the forest canopy on a yearly basis, regardless of their permanence. The distinction between deforestation and forest degradation is made three years after the disturbance occurred by measuring the permanence of the forest disturbance over time. If the forest canopy is disturbed permanently, i.e. shows no signs of forest regrowth over the three years following the disturbance, the ‘forest disturbance’ pixel falls into the deforestation class. If a ‘forest disturbance’ pixel shows clear signs of forest regrowth within the three years following the disturbance, it is classified as forest degradation. In consequence, the distribution of yearly deforestation and forest degradation areas within the measured yearly overall forest disturbance areas are consolidated until 2017, but are estimated (indicated by stars in Figure 21) for the years 2018-2020 on basis of the 16-year average for the period 2002-2017.
The dataset contains yearly information about disturbances in tropical humid forest, starting from year 1984, and about the trajectories over time of each disturbed forest pixel. It differentiates between short-duration and long-duration disturbances, depending on the forest regrowth after the disturbance. Forest regrowth classes only apply to previously deforested areas, not for forest degradation areas. The dataset differentiates between old, young and very young regrowth, depending on when the forest started to grow back after deforestation. Deforestation and forest degradation information are aggregated to the periods ‘before 2000’, ‘2000-2009’ and ‘2010- 2019’ in order to make the map easier to read. However, the data contains the full information, i.e. the precise years of first or consecutive forest disturbances, plus their classification into the deforestation or forest degradation classes
Note: Text from here
from src.processing.jrc_processing import last_of_value, first_of_value
annual_change_maps = xr.concat(objs=[get_jrc_data(lon_min,
lat_min,
lon_max,
lat_max,
dataset="AnnualChange",
year=year
).squeeze().assign_coords({"year": year})
for year in range(1990, 2020)],
dim="year")
%%time
deforestation = rxr.open_rasterio("/gws/nopw/j04/forecol/data/JRC/DeforestationYear/JRC_TMF_DeforestationYear_v1_1982_2019_SAM_ID80_N30_W80.tif").squeeze()
CPU times: user 102 ms, sys: 18.7 ms, total: 121 ms Wall time: 123 ms
annual_change_maps = xr.concat(objs=[rxr.open_rasterio(
f"/gws/nopw/j04/forecol/data/JRC/AnnualChange/tifs/JRC_TMF_AnnualChange_v1_{year}_SAM_ID80_N30_W80.tif"
).squeeze().assign_coords({"year": year})
for year in range(1990, 2020)],
dim="year")
%%time
last_deforestation = xr.zeros_like(deforestation)
CPU times: user 23.2 s, sys: 17 s, total: 40.2 s Wall time: 41.5 s
%%time
last_deforestation.data = last_of_value(annual_change_maps.data, 3)
CPU times: user 9min 51s, sys: 47.1 s, total: 10min 38s Wall time: 10min 38s
%%time
last_deforestation = xr.concat([deforestation.squeeze(), last_deforestation + 1990], dim="observations").max(dim="observations")
CPU times: user 46.2 s, sys: 53.3 s, total: 1min 39s Wall time: 1min 45s
from src.constants import JRC_PATH
%%time
last_deforestation.rio.to_raster(JRC_PATH / "LastDeforestationYear" / "JRC_TMF_LastDeforestationYear_v1_1982_2019_SAM_ID80_N30_W80.tif",
compress='lzw')
CPU times: user 56.1 s, sys: 28.6 s, total: 1min 24s Wall time: 1min 31s
%%time
last_deforestation = xr.zeros_like(deforestation.squeeze())
last_deforestation.data = last_of_value(annual_change_maps.data, 3)
last_deforestation = xr.concat([deforestation.squeeze(), last_deforestation + 1990], dim="observations").max(dim="observations")
<xarray.DataArray (year: 30, y: 37459, x: 37108)>
array([[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]],
[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]],
[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
...
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]],
[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]],
[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]]], dtype=uint8)
Coordinates:
band int64 1
* y (y) float64 30.09 30.09 30.09 30.09 ... 20.0 20.0 20.0 20.0
* x (x) float64 -80.0 -80.0 -80.0 -80.0 ... -70.0 -70.0 -70.0 -70.0
spatial_ref int64 0
* year (year) int64 1990 1991 1992 1993 1994 ... 2016 2017 2018 2019
Attributes:
scale_factor: 1.0
add_offset: 0.0
long_name: Dec1990
grid_mapping: spatial_refarray([[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]],
[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]],
[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
...
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]],
[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]],
[[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
...,
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5],
[5, 5, 5, ..., 5, 5, 5]]], dtype=uint8)array(1)
array([30.094595, 30.094326, 30.094056, ..., 20.000406, 20.000136, 19.999867])
array([-80.000141, -79.999872, -79.999602, ..., -70.000545, -70.000275,
-70.000006])array(0)
array([1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001,
2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013,
2014, 2015, 2016, 2017, 2018, 2019])rxr.open_rasterio("/gws/nopw/j04/forecol/data/JRC/AnnualChange/tifs/JRC_TMF_AnnualChange_v1_2019_SAM_ID80_N30_W80.tif")
<xarray.DataArray (band: 1, y: 37459, x: 37108)>
[1390028572 values with dtype=uint8]
Coordinates:
* band (band) int64 1
* y (y) float64 30.09 30.09 30.09 30.09 ... 20.0 20.0 20.0 20.0
* x (x) float64 -80.0 -80.0 -80.0 -80.0 ... -70.0 -70.0 -70.0 -70.0
spatial_ref int64 0
Attributes:
scale_factor: 1.0
add_offset: 0.0
long_name: Dec2019
grid_mapping: spatial_ref[1390028572 values with dtype=uint8]
array([1])
array([30.094595, 30.094326, 30.094056, ..., 20.000406, 20.000136, 19.999867])
array([-80.000141, -79.999872, -79.999602, ..., -70.000545, -70.000275,
-70.000006])array(0)
i = 0
y, x = regrowth_pixels[i]
regrowth_pixels = np.argwhere(transition_map.data.squeeze() == 30) # 30 = regrowth
# Now I need to get the last deforestation year
first_deforested = np.argwhere(annual_change_maps.data[:, y, x] == 3).min()
print(f"First deforested: {1990 + first_deforested}")
# Now I need to get the last deforestation year
last_deforested = np.argwhere(annual_change_maps.data[:, y, x] == 3).max()
print(f"Last deforested: {1990 + last_deforested}")
First deforested: 1995 Last deforested: 2014
%%time
last_deforestation = xr.zeros_like(deforestation.squeeze())
last_deforestation.data = last_of_value(annual_change_maps.data, 3)
last_deforestation = xr.concat([deforestation.squeeze(), last_deforestation + 1990], dim="observations").max(dim="observations")
CPU times: user 39.1 ms, sys: 61 µs, total: 39.2 ms Wall time: 38.7 ms
np.sum(deforestation > 0)
<xarray.DataArray ()>
array(19750)
Coordinates:
band int64 1
spatial_ref int64 0array(19750)
array(1)
array(0)
np.sum(last_deforestation > 0)
<xarray.DataArray ()>
array(19750)
Coordinates:
band int64 1
spatial_ref int64 0array(19750)
array(1)
array(0)
%%time
first_of_value(annual_change_maps.data, 3)
CPU times: user 32.7 ms, sys: 4.63 ms, total: 37.3 ms Wall time: 37.3 ms
array([[-9999, 26, -9999, ..., 8, 8, 14],
[-9999, 26, -9999, ..., 8, 8, 14],
[-9999, 26, -9999, ..., 14, 14, 14],
...,
[-9999, -9999, -9999, ..., -9999, -9999, -9999],
[-9999, -9999, -9999, ..., -9999, -9999, -9999],
[-9999, -9999, -9999, ..., -9999, -9999, -9999]])
%%time
last_deforested_year = xr.zeros_like(transition_map.squeeze(), dtype=np.uint16)
for y in range(len(annual_change_maps.y)):
for x in range(len(annual_change_maps.x)):
deforested_years = np.argwhere(annual_change_maps.data[:, y, x] == 3)
if len(deforested_years) > 0:
last_deforested_year[y, x] = 1990 + deforested_years.max()
last_deforested_year = xr.concat([last_deforested_year, deforestation.squeeze()], dim="disturbance").max(axis=0)
CPU times: user 3.77 s, sys: 141 ms, total: 3.91 s Wall time: 3.85 s
np.argwhere(annual_change_maps.where(annual_change_maps == 3)[:, y:y+3 , x:x+3].data == 3)
array([[ 5, 0, 0],
[ 5, 0, 1],
[ 5, 0, 2],
[ 5, 1, 2],
[ 5, 2, 0],
[ 5, 2, 1],
[ 5, 2, 2],
[ 6, 0, 0],
[ 6, 0, 1],
[ 6, 0, 2],
[ 6, 1, 2],
[ 6, 2, 0],
[ 6, 2, 1],
[ 6, 2, 2],
[ 7, 0, 0],
[ 7, 0, 1],
[ 7, 0, 2],
[ 7, 1, 2],
[ 7, 2, 0],
[ 7, 2, 1],
[ 7, 2, 2],
[ 8, 0, 0],
[ 8, 0, 1],
[ 8, 0, 2],
[ 8, 1, 2],
[ 8, 2, 0],
[ 8, 2, 1],
[ 8, 2, 2],
[ 9, 0, 0],
[ 9, 0, 1],
[ 9, 0, 2],
[ 9, 1, 2],
[ 9, 2, 0],
[ 9, 2, 1],
[ 9, 2, 2],
[10, 0, 0],
[10, 0, 1],
[10, 0, 2],
[10, 1, 2],
[10, 2, 0],
[10, 2, 1],
[10, 2, 2],
[11, 0, 0],
[11, 0, 1],
[11, 0, 2],
[11, 1, 2],
[11, 2, 0],
[11, 2, 1],
[11, 2, 2],
[12, 0, 0],
[12, 0, 1],
[12, 0, 2],
[12, 1, 2],
[12, 2, 0],
[12, 2, 1],
[12, 2, 2],
[13, 0, 0],
[13, 0, 1],
[13, 0, 2],
[13, 1, 2],
[13, 2, 0],
[13, 2, 1],
[13, 2, 2],
[14, 0, 0],
[14, 0, 1],
[14, 0, 2],
[14, 1, 2],
[14, 2, 0],
[14, 2, 1],
[14, 2, 2],
[15, 0, 0],
[15, 0, 1],
[15, 0, 2],
[15, 1, 2],
[15, 2, 0],
[15, 2, 1],
[15, 2, 2],
[16, 0, 0],
[16, 0, 1],
[16, 0, 2],
[16, 1, 2],
[16, 2, 0],
[16, 2, 1],
[16, 2, 2],
[17, 0, 0],
[17, 0, 1],
[17, 0, 2],
[17, 1, 2],
[17, 2, 0],
[17, 2, 1],
[17, 2, 2],
[18, 0, 0],
[18, 0, 1],
[18, 0, 2],
[18, 1, 2],
[18, 2, 0],
[18, 2, 1],
[18, 2, 2],
[19, 0, 0],
[19, 0, 1],
[19, 0, 2],
[19, 1, 2],
[19, 2, 0],
[19, 2, 1],
[19, 2, 2],
[20, 0, 0],
[20, 0, 1],
[20, 0, 2],
[20, 1, 2],
[20, 2, 0],
[20, 2, 1],
[20, 2, 2],
[21, 0, 0],
[21, 0, 1],
[21, 0, 2],
[21, 1, 2],
[21, 2, 0],
[21, 2, 1],
[21, 2, 2],
[22, 0, 0],
[22, 0, 1],
[22, 0, 2],
[22, 1, 2],
[22, 2, 0],
[22, 2, 1],
[22, 2, 2],
[23, 0, 0],
[23, 0, 1],
[23, 0, 2],
[23, 1, 2],
[23, 2, 0],
[23, 2, 1],
[23, 2, 2],
[24, 0, 0],
[24, 0, 1],
[24, 0, 2],
[24, 1, 2],
[24, 2, 0],
[24, 2, 1],
[24, 2, 2],
[25, 0, 1],
[25, 0, 2],
[25, 1, 2],
[25, 2, 0],
[25, 2, 1],
[25, 2, 2],
[26, 0, 1],
[26, 0, 2],
[26, 1, 2],
[26, 2, 0],
[26, 2, 1],
[26, 2, 2],
[27, 0, 1],
[27, 0, 2],
[27, 1, 2],
[27, 2, 0],
[27, 2, 1],
[27, 2, 2],
[28, 0, 1],
[28, 0, 2],
[28, 1, 2],
[28, 2, 0],
[28, 2, 1],
[28, 2, 2]])
np.argwhere(annual_change_maps.where(annual_change_maps == 3)[:, y:y+10 , x:x+10].data == 3).max(axis=1).shape
(2331,)
import numba
deforestation.squeeze()
<xarray.DataArray (y: 161, x: 462)>
array([[ 0., 2016., 0., ..., 1998., 1998., 2004.],
[ 0., 2016., 0., ..., 1998., 1998., 2004.],
[ 0., 2016., 0., ..., 2004., 2004., 2004.],
...,
[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.]])
Coordinates:
band int64 1
* y (y) float64 -3.279 -3.279 -3.279 -3.28 ... -3.321 -3.322 -3.322
* x (x) float64 -54.94 -54.94 -54.94 ... -54.82 -54.82 -54.82
spatial_ref int64 0
Attributes:
scale_factor: 1.0
add_offset: 0.0
grid_mapping: spatial_refarray([[ 0., 2016., 0., ..., 1998., 1998., 2004.],
[ 0., 2016., 0., ..., 1998., 1998., 2004.],
[ 0., 2016., 0., ..., 2004., 2004., 2004.],
...,
[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.]])array(1)
array([-3.278806, -3.279075, -3.279345, -3.279614, -3.279884, -3.280153,
-3.280423, -3.280692, -3.280962, -3.281231, -3.281501, -3.28177 ,
-3.28204 , -3.282309, -3.282579, -3.282848, -3.283118, -3.283387,
-3.283657, -3.283926, -3.284196, -3.284465, -3.284735, -3.285004,
-3.285274, -3.285543, -3.285813, -3.286082, -3.286352, -3.286621,
-3.286891, -3.28716 , -3.28743 , -3.287699, -3.287969, -3.288238,
-3.288508, -3.288777, -3.289047, -3.289316, -3.289586, -3.289855,
-3.290125, -3.290394, -3.290664, -3.290933, -3.291203, -3.291472,
-3.291742, -3.292011, -3.292281, -3.29255 , -3.29282 , -3.293089,
-3.293359, -3.293628, -3.293898, -3.294167, -3.294437, -3.294706,
-3.294976, -3.295245, -3.295515, -3.295784, -3.296054, -3.296323,
-3.296593, -3.296862, -3.297132, -3.297401, -3.29767 , -3.29794 ,
-3.298209, -3.298479, -3.298748, -3.299018, -3.299287, -3.299557,
-3.299826, -3.300096, -3.300365, -3.300635, -3.300904, -3.301174,
-3.301443, -3.301713, -3.301982, -3.302252, -3.302521, -3.302791,
-3.30306 , -3.30333 , -3.303599, -3.303869, -3.304138, -3.304408,
-3.304677, -3.304947, -3.305216, -3.305486, -3.305755, -3.306025,
-3.306294, -3.306564, -3.306833, -3.307103, -3.307372, -3.307642,
-3.307911, -3.308181, -3.30845 , -3.30872 , -3.308989, -3.309259,
-3.309528, -3.309798, -3.310067, -3.310337, -3.310606, -3.310876,
-3.311145, -3.311415, -3.311684, -3.311954, -3.312223, -3.312493,
-3.312762, -3.313032, -3.313301, -3.313571, -3.31384 , -3.31411 ,
-3.314379, -3.314649, -3.314918, -3.315188, -3.315457, -3.315727,
-3.315996, -3.316266, -3.316535, -3.316805, -3.317074, -3.317344,
-3.317613, -3.317883, -3.318152, -3.318422, -3.318691, -3.318961,
-3.31923 , -3.3195 , -3.319769, -3.320039, -3.320308, -3.320578,
-3.320847, -3.321117, -3.321386, -3.321656, -3.321925])array([-54.942265, -54.941996, -54.941726, ..., -54.818567, -54.818298,
-54.818028])array(0)
last_deforested_year.where(last_deforested_year > 0)
<xarray.DataArray ()>
array(19750)
Coordinates:
band int64 1
spatial_ref int64 0array(19750)
array(1)
array(0)
(last_deforested_year - deforestation.squeeze()).min()
<xarray.DataArray ()>
array(-1989.)
Coordinates:
band int64 1
spatial_ref int64 0array(-1989.)
array(1)
array(0)
np.sum(deforestation > 0)
<xarray.DataArray ()>
array(19750)
Coordinates:
spatial_ref int64 0array(19750)
array(0)
np.sum(last_deforested_year > 0)
<xarray.DataArray ()>
array(19717)
Coordinates:
band int64 1
spatial_ref int64 0array(19717)
array(1)
array(0)
<xarray.DataArray (year: 29)>
array([1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 4, 4, 4, 4], dtype=uint8)
Coordinates:
band int64 1
y float64 -3.279
x float64 -54.93
spatial_ref int64 0
* year (year) int64 1990 1991 1992 1993 1994 ... 2015 2016 2017 2018
Attributes:
scale_factor: 1.0
add_offset: 0.0
long_name: Dec1990
grid_mapping: spatial_refarray([1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 4, 4, 4, 4], dtype=uint8)array(1)
array(-3.27880587)
array(-54.92771263)
array(0)
array([1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001,
2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013,
2014, 2015, 2016, 2017, 2018])deforestation[:, y,x]
<xarray.DataArray (band: 1)>
array([1995.])
Coordinates:
* band (band) int64 1
y float64 -3.279
x float64 -54.93
spatial_ref int64 0
Attributes:
scale_factor: 1.0
add_offset: 0.0
grid_mapping: spatial_refarray([1995.])
array([1])
array(-3.27880587)
array(-54.92771263)
array(0)
degradation[:, y,x]
<xarray.DataArray (band: 1)>
array([0.])
Coordinates:
* band (band) int64 1
y float64 -3.279
x float64 -54.93
spatial_ref int64 0
Attributes:
scale_factor: 1.0
add_offset: 0.0
grid_mapping: spatial_refarray([0.])
array([1])
array(-3.27880587)
array(-54.92771263)
array(0)
annual_change_maps[:, 1, -2]
<xarray.DataArray (year: 29)>
array([1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3], dtype=uint8)
Coordinates:
band int64 1
y float64 -3.279
x float64 -54.82
spatial_ref int64 0
* year (year) int64 1990 1991 1992 1993 1994 ... 2015 2016 2017 2018
Attributes:
scale_factor: 1.0
add_offset: 0.0
long_name: Dec1990
grid_mapping: spatial_refarray([1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3], dtype=uint8)array(1)
array(-3.27907537)
array(-54.81829782)
array(0)
array([1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001,
2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013,
2014, 2015, 2016, 2017, 2018])transition_map[:, 2,1] # 30 = forest regrowth
<xarray.DataArray (band: 1)>
array([30], dtype=uint8)
Coordinates:
* band (band) int64 1
y float64 -3.279
x float64 -54.94
spatial_ref int64 0
Attributes:
scale_factor: 1.0
add_offset: 0.0
long_name: b1
grid_mapping: spatial_refarray([30], dtype=uint8)
array([1])
array(-3.27934486)
array(-54.94199584)
array(0)
deforestation[:, 2, 1]
<xarray.DataArray (band: 1)>
array([2016.])
Coordinates:
* band (band) int64 1
y float64 -3.279
x float64 -54.94
spatial_ref int64 0
Attributes:
scale_factor: 1.0
add_offset: 0.0
grid_mapping: spatial_refarray([2016.])
array([1])
array(-3.27934486)
array(-54.94199584)
array(0)
degradation[:,2,1]
<xarray.DataArray (band: 1)>
array([2008.])
Coordinates:
* band (band) int64 1
y float64 -3.279
x float64 -54.94
spatial_ref int64 0
Attributes:
scale_factor: 1.0
add_offset: 0.0
grid_mapping: spatial_refarray([2008.])
array([1])
array(-3.27934486)
array(-54.94199584)
array(0)
annual_change_maps[:,2, 1]
<xarray.DataArray (year: 29)>
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2,
2, 2, 2, 2, 3, 4, 4], dtype=uint8)
Coordinates:
band int64 1
y float64 -3.279
x float64 -54.94
spatial_ref int64 0
* year (year) int64 1990 1991 1992 1993 1994 ... 2015 2016 2017 2018
Attributes:
scale_factor: 1.0
add_offset: 0.0
long_name: Dec1990
grid_mapping: spatial_refarray([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2,
2, 2, 2, 2, 3, 4, 4], dtype=uint8)array(1)
array(-3.27934486)
array(-54.94199584)
array(0)
array([1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001,
2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013,
2014, 2015, 2016, 2017, 2018])annual_change_maps[:,0,0]
<xarray.DataArray (year: 28)>
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 2, 2], dtype=uint8)
Coordinates:
band int64 1
y float64 -3.279
x float64 -54.94
spatial_ref int64 0
* year (year) int64 1990 1991 1992 1993 1994 ... 2014 2015 2016 2017
Attributes:
scale_factor: 1.0
add_offset: 0.0
long_name: Dec1990
grid_mapping: spatial_refarray([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 2, 2], dtype=uint8)array(1)
array(-3.27880587)
array(-54.94226533)
array(0)
array([1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001,
2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013,
2014, 2015, 2016, 2017])#survey_year = pd.Series(paisagens_master.index, index=paisagens_master.index).str.findall("\w+_(\d{4})\w*").apply(lambda x: int(x[0]))
#paisagens_master["survey_year"] = survey_year
#eba_master["survey_year"] = eba_master.year
#paisagens_master.drop(columns=["geometry_metrics", "geometry_bbox"]).to_csv(PAISAGENSLIDAR_PATH / "paisagens_master_table.csv")
#gpd.GeoDataFrame(paisagens_master, crs=SIRGAS_BRAZIL, geometry="geometry_metrics").to_feather(PAISAGENSLIDAR_PATH / "paisagens_master_table.feather")
#gpd.GeoDataFrame(paisagens_master, crs=SIRGAS_BRAZIL, geometry="geometry_metrics").to_parquet(PAISAGENSLIDAR_PATH / "paisagens_master_table.parquet")
#eba_master.drop(columns=["geometry_metrics", "geometry_bbox", "geometry_metadata"]).to_csv(EBALIDAR_PATH / "eba_master_table.csv")
#gpd.GeoDataFrame(eba_master, crs=SIRGAS_BRAZIL, geometry="geometry_metrics").to_feather(EBALIDAR_PATH / "eba_master_table.feather")
#gpd.GeoDataFrame(eba_master, crs=SIRGAS_BRAZIL, geometry="geometry_metrics").to_parquet(EBALIDAR_PATH / "eba_master_table.parquet")
e_lon_lat = gpd.GeoSeries(eba_master.geometry_bbox,
crs=SIRGAS_BRAZIL).to_crs(WGS84).bounds.rename(columns={"minx": "min_lon",
"miny": "min_lat",
"maxx": "max_lon",
"maxy": "max_lat"})
p_lon_lat = gpd.GeoSeries(paisagens_master.geometry_bbox,
crs=SIRGAS_BRAZIL).to_crs(WGS84).bounds.rename(columns={"minx": "min_lon",
"miny": "min_lat",
"maxx": "max_lon",
"maxy": "max_lat"})